Example of Calibration¶

Import stuff

import warnings

import pandas as pd
import matplotlib.pyplot as plt

from dask.distributed import Client
from dask_jobqueue import SLURMCluster as Cluster

import nevergrad as ng
from mhpc_project.matsch_b2 import CalibrationModel, Variables, Loss, NGO
from geotopy.utils import date_parser, comparison_plot, ProgressBar, ParametersLogger
from geotopy.measures import KGE

Launch Dask cluster using dask-jobqueue

cluster = Cluster()
cluster.adapt(maximum_jobs=128)

client = Client(cluster)

print(cluster.job_script())

#!/usr/bin/env bash

#SBATCH -J dask-worker
#SBATCH -e /dev/shm/dask-worker-%J.err
#SBATCH -o /dev/shm/dask-worker-%J.out
#SBATCH -p regular2
#SBATCH -n 1
#SBATCH --cpus-per-task=4
#SBATCH --mem=2GB
#SBATCH -t 00:20:00
#SBATCH --hint=nomultithread
#SBATCH --hint=compute_bound
export NUMEXPR_NUM_THREADS=1
export NUMEXPR_MAX_THREADS=1
export OMP_NUM_THREADS=1
export TMPDIR=/dev/shm
/home/scampane/scratch/MHPC-project/.env/bin/python -m distributed.cli.dask_worker tcp://10.6.0.4:40470 --nthreads 4 --memory-limit 2.00GB --name name --nanny --death-timeout 60 --local-directory $TMPDIR --lifetime 18m --lifetime-stagger 1m --interface ib0

Load observations, model, loss function and calibration strategy

observations = pd.read_csv('../data/Matsch B2/obs.csv',
                           na_values=['-9999', '-99.99'],
                           usecols=[0, 7],
                           parse_dates=[0],
                           date_parser=date_parser,
                           index_col=0,
                           squeeze=True)
observations.index.rename('datetime', inplace=True)
model = CalibrationModel('../data/Matsch B2/geotop', run_args={'timeout': 120})
variables = Variables('../data/Matsch B2/variables.csv')
measure = KGE(observations)
loss = Loss(model, variables, measure)
calibration = NGO(loss, budget=4096, num_workers=512)

Evaluate the model and plot the risults

simulation = model()
print(f"Before optimization loss is {measure(simulation)}")
comparison_plot(observations,
                simulation,
                desc='Soil moisture content @ 5cm',
                rel=True)
plt.show()

Before optimization loss is 0.7523503598701219

Run calibration

logger = ParametersLogger(loss.massage)

calibration.optimizer.register_callback('tell', logger)

calibration.optimizer.register_callback('tell', ProgressBar(total=calibration.optimizer.budget))

with warnings.catch_warnings():
    
    warnings.simplefilter("ignore")
    recommendation, best_loss = calibration(executor=client)

(256_w,512)-aCMA-ES (mu_w=131.7,w_1=2%) in dimension 21 (seed=<module 'time' (built-in)>, Mon Nov 23 22:06:31 2020)

Plot results and show table of calibrated values

simulation = model(**recommendation)
print(f"After optimization loss is {best_loss}")
comparison_plot(observations,
                simulation,
                desc='Soil moisture content @ 5cm',
                rel=True)
plt.show()

After optimization loss is 0.49338997022896586

pd.DataFrame.from_dict(recommendation, orient='index', columns=['value'])

	value
NVanGenuchten	1.293041
AlphaVanGenuchten	0.000552
ThetaSat	0.578172
ThetaRes	0.066367
NormalHydrConductivity	0.000276
SoilRoughness	168.386071
LSAI	4.317711
CanopyFraction	0.289856
DecayCoeffCanopy	16.953715
RootDepth	286.314931
MinStomatalRes	67.178791
VegReflectVis	0.136903
VegReflNIR	0.295314
VegTransVis	0.017572
VegTransNIR	0.089866
CanDensSurface	0.932729
SoilAlbVisDry	0.063979
SoilAlbNIRDry	0.329262
SoilAlbVisWet	0.215956
SoilAlbNIRWet	0.305602
SoilEmissiv	0.929224

logger.parallel_coordinate_plot()

HiPlot

Loading HiPlot...